#############################################################################
#############################################################################
#Divisive or Descriptive?: How Americans Understand Critical Race Theory
#Safarpour, Lunz Trujillo, Green, Pippert, Lin, and Druckman.
#Journal of Race, Ethnicity, and Politics (2024)
#############################################################################
#############################################################################
rm(list = ls())
#Load packages.
library(survey) 
library(dplyr)
library(margins)
library(car)
library(survey)
library(haven)
library(tidyverse)
library(scales)
library(interplot)
library(coefplot)
library(stargazer)
library(MASS) 
library(srvyr)

d <-read.csv(file = 'CSP_CRT.csv')

#race subsets.
Whites<-subset(d, race.w=="White")
Asian<-subset(d, race.w=="Asian")
Hispanic<-subset(d, race.w=="Hisp")
Black<-subset(d, race.w=="Black")

#Table 1.
##unweighted sample demographics.
round(prop.table(table(d$recage2))*100,0)
round(prop.table(table(d$partylean))*100,0)
round(prop.table(table(d$race.w))*100,0)
round(prop.table(table(d$educ_2))*100,0)
round(prop.table(table(d$gender_re))*100,0)
round(prop.table(table(d$region_rec))*100,0)
#weighted sample demographics.
dat.w <- svydesign(ids = ~ 1,
                   data = d,
                   weights = d$weight.nat) #weight data.
round(addmargins(prop.table(svytable( ~recage2, design = dat.w))) *
        100,0)
round(addmargins(prop.table(svytable( ~partylean, design = dat.w))) *
        100,0)
round(addmargins(prop.table(svytable( ~race.w, design = dat.w))) *
        100,0)
round(addmargins(prop.table(svytable( ~educ_2, design = dat.w))) *
        100,0)
round(addmargins(prop.table(svytable( ~gender_re, design = dat.w))) *
        100,0)
round(addmargins(prop.table(svytable( ~region_rec, design = dat.w))) *
        100,0)

#Figures 2 and 3 (negative affect results).
m3.2<-lm(crt_exp_DV_BinaryNum~crt_treatment*negativeAffect, data=d)
summary(m3.2)
summary(margins(m3.2, variables="crt_treatment", at=list(negativeAffect=c(0, 0.25, 0.5, 0.75, 1)),type="response"))

#Figure 2.
pdf("Safarpour_etal_Fig2.pdf")
interplot(m=m3.2, var1="crt_treatment", var2="negativeAffect")+
  xlab("Negative Affect Toward Black People")+
  ylab("Estimated Coefficient for 'CRT' Condition\n (VS 'Defintion' Condition)")+
  theme_bw()
dev.off()

d$negativeAffectQuant<-gtools::quantcut(d$negativeAffect, q=4, na.rm=TRUE)
def1<-as.data.frame(addmargins(prop.table(table(d$negativeAffectQuant, d$defsupport),1))*100)
crt1<-as.data.frame(addmargins(prop.table(table(d$negativeAffectQuant, d$crtsupport),1))*100)
def1<-def1[which(def1$Var1!="Sum"),]
def1<-def1[which(def1$Var2!="Sum"),]
crt1<-crt1[which(crt1$Var1!="Sum"),]
crt1<-crt1[which(crt1$Var2!="Sum"),]
names(crt1)[names(crt1) == 'Freq'] <- "CRT"
sum<-cbind(def1, crt1)
sum<-subset(sum, select=-c(1,2))
sum$diff<-sum$CRT-sum$Freq
sum$levelAffect<-rep(c(1,2,3,4),3)
sum$levelAffect<-rep(c("Low (0-0.48)", "Medium-Low (0.49-0.5)",
                       "Medium-High (0.51-0.515)", "High (0.52-1)"),3)
sum$levelAffect2<-factor(sum$levelAffect, ordered=TRUE, levels=c("Low (0-0.48)", "Medium-Low (0.49-0.5)",
                                                                 "Medium-High (0.51-0.515)", "High (0.52-1)"))
sum$Var2<-as.character(sum$Var2)
sum$Var2[sum$Var2=="Neither support not oppose"]<-"Neither support nor oppose"
sum$levelSupport<-factor(sum$Var2, ordered=TRUE, levels=c("Support","Neither support nor oppose", 
                                                          "Oppose"))


pdf("Safarpour_etal_Fig3.pdf")
p<-ggplot(sum, aes(fill=levelAffect2,x=levelSupport, y=diff))+
  geom_bar(position="dodge",stat="identity")+scale_fill_grey()
p+theme_bw()+
  labs(x="Attitude Toward Teaching \nCRT/ Legacy of Racism in Public Schools", y="Percentage Point Difference in Attitude \nToward Teaching CRT VS Legacy of Racism",
       fill="Negative Affect")
dev.off()

#appendix fig b1.
m3.2O<-lm(crt_exp_DV_ord~crt_treatment*negativeAffect, data=d)
summary(m3.2O)
summary(margins(m3.2O, variables="crt_treatment", at=list(negativeAffect=c(0, 0.25, 0.5, 0.75, 1)),type="response"))
pdf("CRTEffectsbyNegativeAffectOrdDV.pdf")
interplot(m=m3.2O, var1="crt_treatment", var2="negativeAffect")+
  xlab("Negative Affect Toward Black People")+
  ylab("Estimated Coefficient for 'CRT' Condition\n (VS 'Defintion' Condition)")+
  theme_bw()
dev.off()
#appendix table b3.
m3.1<-lm(crt_exp_DV_BinaryNum~crt_treatment*partylean, data=d)
summary(m3.1)
m3.1w<-lm(crt_exp_DV_BinaryNum~crt_treatment*partylean, data=d, weights=weight.nat)
summary(m3.1w)
m3.2<-lm(crt_exp_OpposeBinaryNum~crt_treatment*negativeAffect, data=d)
summary(m3.2)
m3.2w<-lm(crt_exp_OpposeBinaryNum~crt_treatment*negativeAffect, data=d, weights=weight.nat)
summary(m3.2w)
m3.20<-lm(crt_exp_DV_ord~crt_treatment*negativeAffect, data=d)
summary(m3.20)
m3.20w<-lm(crt_exp_DV_ord~crt_treatment*negativeAffect, data=d, weights=weight.nat)
summary(m3.20w)
stargazer(m3.1,m3.1w,m3.2,m3.2w,m3.20,m3.20w,
          type = "html", out="CRT Support CATE.html",
          #title="", 
          covariate.labels = c("CRT Condition", "Independent", "Republican/Lean Republican", 
                               "CRT Condition*Independent", "CRT Condition*Republican",
                               "Negative Affect", "CRT Condition*Negative Affect"),
          dep.var.labels=c("Binary Support", "Binary Support", "Ordinal Support (5 pt)"),
          column.labels = c("Unweighted", "Weighted","Unweighted", "Weighted","Unweighted", "Weighted"),
          notes="OLS Estimates. Standard errors in parentheses. Binary dependent variable is coded 1=Strongly/ somewhat support, 0=otherwise. Ordinal dependent variable is coded 0=Strongly oppose, 0.25=Somewhat oppose, 0.5=Neither support nor oppose, 0.75= Somewhat support, 1= Strongly support. CRT condition is coded 1 if respondents were asked `Do you support or oppose teaching Critical Race Theory (CRT) in public schools?' CRT condition is coded as 0 if respondents were asked ' Do you support or oppose teaching how racism continues to impact American society today in public schools?' Negative affect is coded 0 to 1.",
          notes.align = "l", notes.label = "Notes:" )
#Figure 3 (pid results).
results<-as.data.frame(summary(margins(m3.1, variables="crt_treatment", at=list(partylean=c("Democrat/ Lean Democrat", "Independent","Republican/ Lean Republican")),type="response")))
#output: 1=Dems; 2=Independents; 3=Republicans.
options(digits=4)
pdf("Safarpour_etal_Fig4.pdf")
TA<-ggplot(d=results, aes(partylean, AME))+
  geom_hline(yintercept=0, linetype="dashed", color="gray28")+
  geom_point(position = "dodge",
             stat = "identity")+
  geom_errorbar(aes(ymin=lower, ymax=upper), lty="twodash", width=0.7)+
  labs(y="Estimated Coefficient for 'CRT' Condition\n (VS 'Defintion' Condition)", 
       x="Party Identification"
  )+  theme_bw()+
  geom_text(
    aes(label = round((AME), 2)),
    position = position_dodge(width = 1.5),
    vjust = -0.5,
    hjust=-0.25, #shift text horizontally to the right.
    size = 3)
TA
dev.off()

#toplines reported in text.
round(addmargins(prop.table(table(d$defsupport, d$partylean),2))*100,0)
round(addmargins(prop.table(table(d$crtsupport, d$partylean),2))*100,0)
d$crt_exp_DVO<-as.factor(d$crt_exp_DVO)
oprobit<-polr(crt_exp_DVO ~ crt_treatment*partylean, data=d, method='probit', Hess=TRUE)

#ordered results for appendix.
sample <- d[complete.cases(d$crt_exp_DVO, d$crt_treatment, d$partylean)==T,]
preds <- predict(oprobit, sample, type='probs')
summary(preds) #baseline predicted probabilities, one for each cut point.
newdat <- cbind(sample, predict(oprobit, sample, type = "probs"))
newdat <- subset(newdat, select=c("crt_exp_DVO", "crt_treatment", "partylean", "Oppose", "Neither support not oppose", "Support"))
library(reshape2)
lnewdat <- melt(newdat, id.vars = c("crt_exp_DVO", "crt_treatment", "partylean"),
                variable.name = "Level", value.name="Probability")
lnewdat$`Party Identification`<-lnewdat$partylean
lnewdat$Condition<-"Definition"
lnewdat$Condition[lnewdat$crt_treatment==1]<-"CRT"
lnewdat$Level2<-recode_factor(lnewdat$Level, Oppose="Oppose",`Neither support not oppose`= "Neither support\n nor oppose", Support="Support")
#appendix figure b2.
png("OProbit.png")
ggplot(lnewdat, aes(x = Condition, y = Probability, colour = `Party Identification`, shape=`Party Identification`)) +
  geom_point() +  facet_wrap(~Level2, nrow=1) +theme_bw()+labs(y="Predicted Probability")
dev.off()

#appendix table b2.
m1<-lm(crt_exp_DV_BinaryNum~crt_treatment*race.w, data=d)
m1.w<-lm(crt_exp_DV_BinaryNum~crt_treatment*race.w, data=d, weights=weight.nat)
m11<-lm(crt_exp_DV_BinaryNum~crt_treatment*partylean, data=Whites)
m11.w<-lm(crt_exp_DV_BinaryNum~crt_treatment*partylean, data=Whites, weights=weight.nat)

m2<-lm(crt_exp_DV_BinaryNum~crt_treatment*partylean, data=Black)
m2.w<-lm(crt_exp_DV_BinaryNum~crt_treatment*partylean, data=Black, weights=weight.nat)

m3<-lm(crt_exp_DV_BinaryNum~crt_treatment*partylean, data=Hispanic)
m3.w<-lm(crt_exp_DV_BinaryNum~crt_treatment*partylean, data=Hispanic, weights=weight.nat)

m4<-lm(crt_exp_DV_BinaryNum~crt_treatment*partylean, data=Asian)
m4.w<-lm(crt_exp_DV_BinaryNum~crt_treatment*partylean, data=Asian, weights=weight.nat)

stargazer(m1, m1.w, m11, m11.w, m2, m2.w, m3, m3.w, m4, m4.w,
          type = "html", out="CRT Support Race.html",
          covariate.labels = c("CRT Condition", "Black", "Hispanic", "Other Race", "White",
                               "CRT Condition*Black", "CRT Condition*Hispanic", 
                               "CRT Condition*Other", "CRT Condition*White",
                               "Independent", "Republican/ Lean Republican", 
                               "CRT Condition*Independent", "CRT Condition*Republican"),
          dep.var.labels=c( "Binary Support"),
          column.labels = c("Unweighted", "Weighted", "Whites: Unweighted", "Whites: Weighted", 
                            "Blacks: Unweighted", "Blacks: Weighted",
                            "Hisp.: Unweighted", "Hisp.: Weighted", "Asian: Unweighted", "Asian: Weighted"),
          notes="OLS Estimates. Standard errors in parentheses. Binary dependent variable is coded 1=Strongly/ somewhat support, 0=otherwise. CRT condition is coded 1 if respondents were asked `Do you support or oppose teaching Critical Race Theory (CRT) in public schools?' CRT condition is coded as 0 if respondents were asked ' Do you support or oppose teaching how racism continues to impact American society today in public schools?'",
          notes.align = "l", notes.label = "Notes:" , single.row = F)
# Calculate marginal effects
m1<-lm(crt_exp_DV_BinaryNum~crt_treatment*race.w, data=d)
margins_table <- summary(margins(m1, variables = c("crt_treatment"), at = list(race.w = c("White", "Black", "Hisp", "Asian", "Other"))))
margins_table$Race<-c("White", "Black", "Hispanic", "Asian", "Other")

# Specify the order of levels for the x-axis
margins_table$Race <- factor(margins_table$Race, levels = c("White", "Black", "Hispanic", "Asian", "Other"))

margins_table$partylean<-"Overall"
margins_table <- margins_table %>%
  select(-race.w)
tab1 <- summary(margins(m11, variables = c("crt_treatment"), at = list(partylean = c("Democrat/ Lean Democrat", "Independent", "Republican/ Lean Republican"))))
tab1$Race<-"White"
tab2 <- summary(margins(m2, variables = c("crt_treatment"), at = list(partylean = c("Democrat/ Lean Democrat", "Independent", "Republican/ Lean Republican"))))
tab2$Race<-"Black"
tab3 <- summary(margins(m3, variables = c("crt_treatment"), at = list(partylean = c("Democrat/ Lean Democrat", "Independent", "Republican/ Lean Republican"))))
tab3$Race<-"Hispanic"
tab4 <- summary(margins(m4, variables = c("crt_treatment"), at = list(partylean = c("Democrat/ Lean Democrat", "Independent", "Republican/ Lean Republican"))))
tab4$Race<-"Asian"
sum<-rbind(margins_table,tab1,tab2,tab3,tab4)

#figure 5.
pdf("Safarpour_etal_Fig5.pdf")
ggplot(sum, aes(x=Race , y = AME)) +
  facet_wrap(~partylean, nrow=2)+
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
  labs(x = "Party Identification", y = "Estimated Coefficient for CRT Condition") +
  theme_bw() +
  scale_y_continuous(labels = scales::number_format(scale = 1, accuracy = 0.01))+
  geom_hline(yintercept = 0, linetype = "dashed")
dev.off()

#Figure 6.
`Pre-Registered Model`<-lm(crt_familiarBinaryNum~partylean+ideol3+nonwhite+age+gender_re+parent+`educ_2_High.School.or.Less`, data=d) 
`Robustness Check`<-lm(crt_familiarBinaryNum~partylean+ideol3+nonwhite+age+gender_re+parent+`educ_2_High.School.or.Less`+crt_treatment, data=d) 

pdf("Safarpour_etal_Fig6.pdf")
coefplot::multiplot(`Pre-Registered Model`, `Robustness Check`)+ 
  ggtitle("")+xlab("Coefficient Value")+
  theme_bw()+
  scale_y_discrete(labels=c("Intercept","Independent", "Republican/Lean Republican", 
                            "Ideology: Moderate", "Ideology: Conservative", 
                            "Nonwhite", "Age", "Male",  
                            "Parent", "High School or Less", "CRT Treatment"))
dev.off()

#Appendix figure b3 and appendix table b4.
m4<-lm(crt_exp_DV_BinaryNum~crt_treatment+negativeAffect+partylean+ideol3+Trumpvoter+nonwhite+gender_re+antifaJan6Binary+migrantsCovidMythBinary+parent+Vaccinated, data=d) #need to add  and belief in covid-19 myths to predictors.
m4.w<-lm(crt_exp_DV_BinaryNum~crt_treatment+negativeAffect+partylean+ideol3+Trumpvoter+nonwhite+gender_re+antifaJan6Binary+migrantsCovidMythBinary+parent+Vaccinated, data=d, weights=weight.nat) #need to add  and belief in covid-19 myths to predictors.

stargazer(m4,m4.w,
          type = "html", out="CRT Support Predictors.html",
          covariate.labels = c("CRT Condition", "Negative Affect", "Independent", "Republican/Lean Republican", 
                               "Ideology: Moderate", "Ideology: Conservative", "Trump supporter",
                               "Nonwhite", "Male", "Jan 6 was Antifa", "Migrants Cause COVID Surge",
                               "Parent", "Vaccinated for COVID-19"),
          dep.var.labels=c( "Support Teaching CRT/ Legacy of Racism In Public Schools"),
          column.labels = c("Unweighted", "Weighted"),
          notes="OLS Estimates. Standard errors in parentheses. Binary dependent variable is coded 1=Strongly/ somewhat support, 0=otherwise. CRT condition is coded 1 if respondents were asked `Do you support or oppose teaching Critical Race Theory (CRT) in public schools?' CRT condition is coded as 0 if respondents were asked ' Do you support or oppose teaching how racism continues to impact American society today in public schools?' Negative affect is coded 0 to 1. Antifa and Migrants are coded 0=Innacurate, 1=Accurate/ Not sure. Election Stolen coded 0=disagree/neither, 1= agree. All independent variables are coded between 0 and 1.",
          notes.align = "l", notes.label = "Notes:" , single.row = T)

m4<-lm(crt_exp_DV_BinaryNum~crt_treatment+negativeAffect+partylean+ideol3+Trumpvoter+nonwhite+gender_re+antifaJan6Binary+migrantsCovidMythBinary+parent+Vaccinated, data=d) #need to add  and belief in covid-19 myths to predictors.
summary(m4)
png("PredictorsSupportTeachingCRT.png")
coefplot::coefplot(m4)+ 
  ggtitle("Support for Teaching CRT/Legacy of Racism in Schools")+
  xlab("Coefficient Value")+ theme_bw()+
  scale_y_discrete(labels=c("Intercept","CRT Treatment","Negative Affect","Independent", "Republican/Lean Republican", 
                            "Ideology: Moderate", "Ideology: Conservative", 
                            "Trump supporter","Nonwhite",  "Male",  
                            "Jan 6 was Antifa", "Migrants Cause COVID Surge",
                            "Parent", "Vaccinated for COVID-19"))
dev.off()

#appendix table b5.
#predictors of how well CRT describes American society--
m5<-lm(crt_societyBinaryNum~partylean+ideol3+nonwhite, data=d) 
m5w<-lm(crt_societyBinaryNum~partylean+ideol3+nonwhite, data=d, weights=weight.nat) 
m5.1<-lm(crt_societyBinaryNum~partylean+ideol3+nonwhite+age+gender_re+negativeAffect+Trumpvoter+antifaJan6Binary+migrantsCovidMythBinary+parent+Vaccinated, data=d) 
m5.1w<-lm(crt_societyBinaryNum~partylean+ideol3+nonwhite+age+gender_re+negativeAffect+
            Trumpvoter+antifaJan6Binary+migrantsCovidMythBinary+parent+Vaccinated, data=d, weights=weight.nat) 
stargazer(m5,m5w,m5.1, m5.1w,
          type = "html", out="CRT Society Predictors.html",
          covariate.labels = c(  "Independent", "Republican/Lean Republican", 
                                 "Ideology: Moderate", "Ideology: Conservative", 
                                 "Nonwhite", "Age", "Male", "Negative Affect", "Trump supporter",
                                 "Jan 6 was Antifa", "Migrants Cause COVID Surge",
                                 "Parent", "Vaccinated for COVID-19"),
          dep.var.labels=c( "CRT Describes Society Extremely/Very Well"),
          column.labels = c("Unweighted", "Weighted", "Unweighted", "Weighted"),
          notes="OLS Estimates. Standard errors in parentheses. Binary dependent variable is coded 1=CRT describes society well, 0=otherwise. Negative affect is coded 0 to 1. Antifa and Migrants are coded 0=Innacurate, 1=Accurate/ Not sure. All independent variables are coded between 0 and 1.",
          notes.align = "l", notes.label = "Notes:" , single.row = T)

#figure 7.
code1<-as.data.frame(round((prop.table(table(d$Code1R))*100),1))
code2<-as.data.frame(round((prop.table(table(d$Code2R))*100),1))
code3<-as.data.frame(round((prop.table(table(d$Code3R))*100),1))

code1$Code<-"Code 1"
code2$Code<-"Code 2"
code3$Code<-"Code 3"

sum<-rbind(code1, code2, code3)
sum$Var2<-c("College", "Don't\nKnow",  "Kids", "Negative", "Neutral/\nPositive","No Teach", "Other", "Partisan", "Wikipedia",
            "College", "Don't\nKnow",  "Kids", "Negative", "Neutral/\nPositive","No Teach", "Other", "Partisan", "Wikipedia",
            "College", "Don't\nKnow",  "Kids", "Negative", "Neutral/\nPositive","No Teach", "Other", "Partisan")

pdf("Safarpour_etal_Fig7.pdf")
ggplot(sum, aes(x=Var2, y=Freq))+
  geom_bar(stat="identity")+
  facet_wrap(~Code, ncol=1)+
  theme_bw()+theme(axis.text.x = element_text(angle = 0))+
  labs(x="In your own words, please describe what Critical Race Theory (CRT) means to you", y="Percent" )+
  geom_text(
    aes(label = Freq),
    position = position_dodge(width = 0),
    vjust = -0.4,
    hjust=0.4, #shift text horizontally to the right.
    size = 3
  )+ylim(0,65)
dev.off()

#crt meaning by experimental condition.
#regressions to confirm small/null differences (reported in text)
m1<-lm(Negative~crt_treatment, data=d)
summary(m1)
m2<-lm(NeutralPos~crt_treatment, data=d)
summary(m2)
m3<-lm(DK~crt_treatment, data=d)
summary(m3)

#racism/ black open ended tretament.
d$BlkTreatRacism<-ifelse((d$Code1=="RACISM"| d$Code1=="BLK/MINORITY TREAT"),1,0)
table(d$BlkTreatRacism, d$Code1)
d$BlkTreatRacism[d$Code2=="RACISM"| d$Code2=="BLK/MINORITY TREAT"]<-1
d$BlkTreatRacism[d$Code3=="RACISM"| d$Code3=="BLK/MINORITY TREAT"]<-1
m4<-lm(BlkTreatRacism~crt_treatment, data=d)
summary(m4)

#Weighted Crosstabs.
d.w <- svydesign(ids = ~ 1,
                 data = d,
                 weights = d$weight.nat) #weight data.

code1<-as.data.frame(round((prop.table(svytable( ~Code1R, design = d.w))*100),1))
code2<-as.data.frame(round((prop.table(svytable( ~Code2R, design = d.w))*100),1))
code3<-as.data.frame(round((prop.table(svytable( ~Code3R, design = d.w))*100),1))
code1$Code<-"Code 1"
code2$Code<-"Code 2"
code3$Code<-"Code 3"
code1<-code1%>% rename(Var1=Code1R)
code2<-code2%>% rename(Var1=Code2R)
code3<-code3%>% rename(Var1=Code3R)
sum<-rbind(code1, code2, code3)

#appendix figure b4.
png("WeightedCRTCodes.png")
ggplot(sum, aes(x=Var1, y=Freq))+
  geom_bar(stat="identity")+
  facet_wrap(~Code, ncol=1)+
  theme_bw()+theme(axis.text.x = element_text(angle = 90))+
  labs(x="Code", y="Percent" )+ggtitle("Weighted, Collapsed CRT Meaning Codes")+
  geom_text(
    aes(label = Freq),
    position = position_dodge(width = 0),
    vjust = -0.4,
    hjust=0.4, #shift text horizontally to the right.
    size = 3
  )+ylim(0,70)
dev.off()

#CRT meaning by partisanship.-----
m1<-lm(DVSCRT~partylean, data=d) 
preds1<-as.data.frame(summary(prediction(m1, calculate_se=TRUE, at = list(partylean = c("Republican/ Lean Republican", "Democrat/ Lean Democrat")))))
preds1$Code<-"Divisive"
rm(m1)
m1<-lm(BlkVictimCRT~partylean, data=d) 
preds2<-as.data.frame(summary(prediction(m1, calculate_se=TRUE, at = list(partylean = c("Republican/ Lean Republican", "Democrat/ Lean Democrat")))))
preds2$Code<-"Black Victim/ Favors"
rm(m1)
m1<-lm(KidsCRT~partylean, data=d) 
preds3<-as.data.frame(summary(prediction(m1, calculate_se=TRUE, at = list(partylean = c("Republican/ Lean Republican", "Democrat/ Lean Democrat")))))
preds3$Code<-"Kids"
rm(m1)
m1<-lm(PartisanCRT~partylean, data=d) 
preds4<-as.data.frame(summary(prediction(m1, calculate_se=TRUE, at = list(partylean = c("Republican/ Lean Republican", "Democrat/ Lean Democrat")))))
preds4$Code<-"Partisan"
rm(m1)
m1<-lm(RacistAWCRT~partylean, data=d) 
preds5<-as.data.frame(summary(prediction(m1, calculate_se=TRUE, at = list(partylean = c("Republican/ Lean Republican", "Democrat/ Lean Democrat")))))
preds5$Code<-"Racist or Anti-White"

rm(m1)
m1<-lm(PROPCRT~partylean, data=d) 
preds6<-as.data.frame(summary(prediction(m1, calculate_se=TRUE, at = list(partylean = c("Republican/ Lean Republican", "Democrat/ Lean Democrat")))))
preds6$Code<-"Propaganda"
rm(m1)
preds<-rbind(preds1, preds2, preds3, preds4, preds5, preds6)
preds$Prediction<-preds$Prediction*100
preds$lower<-preds$lower*100
preds$upper<-preds$upper*100
preds$Party<-rep(c("Republican/ \nLean Republican","Democrat/ \nLean Democrat"),6)

#Figure 8.
pdf("Safarpour_etal_Fig8.pdf")
p<-ggplot(preds, aes( y=Prediction, x=Party,  fill=Party)) + facet_wrap(~Code)+
  geom_bar(position = position_dodge(width = 0.8), #set doge wider than bar width to create space between bars.
           stat = "identity",
           width = 0.7) + theme_bw()+
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.2,
                position=position_dodge(0.8)) +
  labs(x="Party Identification", 
       y="Percent")+
  theme(legend.position = "none")+scale_fill_grey()
p+scale_y_continuous( breaks=c(0,5,10,15,20),labels=c("0%","5%" ,"10%","15%","20%"))
dev.off() 

#negative, neutral/positive, dk, etc. by party.
prop.table(table(d$Negative, d$partylean),2)*100
prop.table(table(d$DK, d$partylean),2)*100
prop.table(table(d$NeutralPos, d$partylean),2)*100


m1<-lm(Negative~partylean, data=d) 
summary(m1)
preds1<-as.data.frame(summary(prediction(m1, calculate_se=TRUE, at = list(partylean = c("Republican/ Lean Republican", "Democrat/ Lean Democrat")))))
preds1$Code<-"Negative"
rm(m1)
m1<-lm(NeutralPos~partylean, data=d) 
summary(m1)
preds2<-as.data.frame(summary(prediction(m1, calculate_se=TRUE, at = list(partylean = c("Republican/ Lean Republican", "Democrat/ Lean Democrat")))))
preds2$Code<-"Neutral/Positive"
rm(m1)
m1<-lm(DK~partylean, data=d) 
summary(m1)
preds3<-as.data.frame(summary(prediction(m1, calculate_se=TRUE, at = list(partylean = c("Republican/ Lean Republican", "Democrat/ Lean Democrat")))))
preds3$Code<-"Don't Know"
preds<-rbind(preds1, preds2, preds3)
preds$Prediction<-preds$Prediction*100
preds$lower<-preds$lower*100
preds$upper<-preds$upper*100
preds$Party<-rep(c("Republican/ \nLean Republican","Democrat/ \nLean Democrat"),3)
p<-ggplot(preds, aes( y=Prediction, x=Party,  fill=Party)) + facet_wrap(~Code)+
  geom_bar(position = position_dodge(width = 0.8), 
           stat = "identity", width = 0.7) + theme_bw()+
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.2,
                position=position_dodge(0.8)) +
  labs(x="Party Identification",  y="Percent")+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("mediumblue",  "red"))
p+scale_y_continuous( breaks=c(0,10,20,30,40,50),labels=c("0%" ,"10%","20%", "30%", "40%", "50%"))

#Appendix OE analysis by race.
code1<-as.data.frame(round((prop.table(table( d$Code1R, d$race.w),2)*100),1))
code2<-as.data.frame(round((prop.table(table( d$Code2R, d$race.w),2)*100),1))
code3<-as.data.frame(round((prop.table(table( d$Code3R, d$race.w),2)*100),1))
code1$Code<-"Code 1"
code2$Code<-"Code 2"
code3$Code<-"Code 3"
sum<-rbind(code1, code2, code3)
sum <- sum[sum$Var2 != "Other", ]
names(sum)[names(sum) == "Var2"] <- "Race"
#appendix figure b5.
png("crtcodesrace.png")
ggplot(sum, aes(x=Var1, y=Freq, fill=Race))+
  geom_bar(stat = "identity", position = position_dodge(width = 0.7), width = 0.6) +
  facet_wrap(~Code, ncol=1)+
  theme_bw()+theme(axis.text.x = element_text(angle = 90, size=8))+
  labs(x="Code", y="Percent" )+
  geom_text(
    aes(label = round(Freq,0)),
    position = position_dodge(width = 0.7),
    vjust = -0.4,
    hjust=0.4, #shift text horizontally to the right.
    size = 2
  )+ scale_fill_grey()
dev.off()

#Fig 9. White feeling therm by party/ anti-white/divisive open-ended responses.
m1<-lm(therm1_6~RacistAWCRT, data=subset(d, race.w=="White"));summary(m1)
preds3<-as.data.frame(summary(prediction(m1, calculate_se=TRUE, at = 
                                           list(RacistAWCRT=c(0,1)))))
preds3$OE<-c("No", "Yes")
pdf("Safarpour_etal_Fig9.pdf")
ggplot(preds3, aes( y=Prediction, x=OE)) + 
  geom_point(position = "dodge", stat = "identity")+
  geom_errorbar(aes(ymin=lower, ymax=upper), lty="twodash", width=0.3)+
  labs(x="Defined CRT as Racist or Anti-White\n in Open-Ended Response", 
       y="Average White Feeling Thermometer")+theme_bw()
dev.off() 

#Overall Opinions.
m1<-lm(crt_exp_DV_ord~crt_treatment, data=d)
m1.1<-lm(crt_exp_DV_ord~crt_treatment, data=d, weights= weight.nat)
m2<-lm(crt_exp_DV_BinaryNum~crt_treatment, data=d)
m2.1<-lm(crt_exp_DV_BinaryNum~crt_treatment, data=d, weights= weight.nat)
#appendix table b1. 
stargazer(m1,m1.1,m2,m2.1,
          type = "html", out="CRT Support Basic.html",
          covariate.labels = c("CRT Condition"),
          dep.var.labels=c( "Ordinal Support (5 pt)", "Binary Support"),
          column.labels = c("Unweighted", "Weighted", "Unweighted", "Weighted"),
          notes="OLS Estimates. Standard errors in parentheses. Ordinal dependent variable is coded 0=Strongly oppose, 0.25=Somewhat oppose, 0.5=Neither support nor oppose, 0.75= Somewhat support, 1= Strongly support. Binary dependent variable is coded 1=Strongly/ somewhat support, 0=otherwise. CRT condition is coded 1 if respondents were asked `Do you support or oppose teaching Critical Race Theory (CRT) in public schools?' CRT condition is coded as 0 if respondents were asked ' Do you support or oppose teaching how racism continues to impact American society today in public schools?'",
          notes.align = "l", notes.label = "Notes:" , single.row = F)

#figure 10.
# Group the data by race.w and calculate average crt_exp_DV_BinaryNum for each crt_treatment value
grouped_data <- d %>%
  group_by(race.w, crt_treatment) %>%
  summarize(Avg_crt_exp_DV_BinaryNum = mean(crt_exp_DV_BinaryNum, na.rm = TRUE))
# Print the grouped data
print(grouped_data)
# Recode crt_treatment values and rename columns
grouped_data <- grouped_data %>%
  mutate(Condition = ifelse(crt_treatment == 0, "Definition", "CRT")) %>%
  select(race.w, Condition, Avg_crt_exp_DV_BinaryNum)
# Calculate the overall average for each value of crt_condition
overall_avg <- d %>%
  group_by(crt_treatment) %>%
  summarize(Avg_crt_exp_DV_BinaryNum = mean(crt_exp_DV_BinaryNum, na.rm = TRUE)) %>%
  mutate(race.w = "Overall")
overall_avg <-overall_avg %>%
  mutate(Condition = ifelse(crt_treatment == 0, "Definition", "CRT")) %>%
  select(race.w, Condition, Avg_crt_exp_DV_BinaryNum)
# Combine the dataset with overall_avg
grouped_data <- bind_rows(grouped_data, overall_avg)
# Specify the order of levels for the x-axis
grouped_data$race.w <- factor(grouped_data$race.w, levels = c("Overall","White", "Black", "Hisp", "Asian", "Other"))
grouped_data$Avg_crt_exp_DV_BinaryNum<-grouped_data$Avg_crt_exp_DV_BinaryNum*100
#bar plot.
pdf("Safarpour_etal_Fig10.pdf")
ggplot(grouped_data, aes(x = race.w, y = Avg_crt_exp_DV_BinaryNum, fill = Condition)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.7), width = 0.6) +
  geom_text(aes(label = paste0(round(Avg_crt_exp_DV_BinaryNum, 0), "%"), y = Avg_crt_exp_DV_BinaryNum + 0.01),
            position = position_dodge(width = 0.7), vjust = -0.5) +
  labs(x = "Race/Ethnicity", y = "Percent Support") +
  theme_bw() +
  scale_fill_manual(values = c("Definition" = "gray", "CRT" = "black")) +
  guides(fill = guide_legend(title = "Condition"))+
  scale_x_discrete(labels = c("Overall","White", "Black", "Hispanic", "Asian", "Other"))+
  scale_y_continuous(labels = scales::percent_format(scale = 1, accuracy = 1))
dev.off()

#RACE interactive model (overall).
m1<-lm(crt_exp_DV_BinaryNum~crt_treatment*race.w, data=d);summary(m1)
# Calculate marginal effects
margins_table <- summary(margins(m1, variables = c("crt_treatment"), at = list(race.w = c("White", "Black", "Hisp", "Asian", "Other"))))
margins_table$Race<-c("White", "Black", "Hispanic", "Asian", "Other")
# Specify the order of levels for the x-axis
margins_table$Race <- factor(margins_table$Race, levels = c("White", "Black", "Hispanic", "Asian", "Other"))
# Plot.
ggplot(margins_table, aes(x = Race, y = AME)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
  labs(x = "Race/Ethnicity", y = "Estimated Coefficient for CRT Condition") +
  theme_bw() +
  scale_y_continuous(labels = scales::number_format(scale = 1, accuracy = 0.01))

#endnote xii- fox news analysis.
famPredict<-lm(crt_familiarBinaryNum~partylean+ideol3+nonwhite+age+gender_re+parent+`educ_2_High.School.or.Less`+FoxViewer, data=d) 
summary(famPredict)
negativeOE<-lm(Negative~partylean+ideol3+nonwhite+age+gender_re+parent+`educ_2_High.School.or.Less`+FoxViewer, data=d) 
summary(negativeOE)


#Table 2 and Table 3 (Breakdowns of larger categories).
addmargins(round((prop.table(table( d$Code1, d$Code1R),2)*100),0))


#Read in ICR data (note: data is subset to just those that were coded by two coders- responses were randomly selected for dual coders.)
rm(list = ls())

sum<-readxl::read_excel("ICRdata.xlsx", sheet="Sheet1", 
                        col_names=TRUE)
#Calculate ICR.
sum$AgreeCode1<-0
sum$AgreeCode1[sum$Code1==sum$Code1x|
                 sum$Code1==sum$Code2x|
                 sum$Code1==sum$Code3x]<-1

sum$AgreeCode2<-0
sum$AgreeCode2[sum$Code2==sum$Code1x|
                 sum$Code2==sum$Code2x|
                 sum$Code2==sum$Code3x]<-1
sum$AgreeCode2[is.na(sum$Code2) & is.na(sum$Code2x)]<-NA
table(sum$AgreeCode2, exclude=NULL)
table(sum$Code2, exclude=NULL)
table(sum$Code2x, exclude=NULL)

sum$AgreeCode3<-0
sum$AgreeCode3[sum$Code3==sum$Code1x|
                 sum$Code3==sum$Code2x|
                 sum$Code3==sum$Code3x]<-1
sum$AgreeCode3[is.na(sum$Code3) & is.na(sum$Code3x)]<-NA

sum$AgreeScore<- 0
sum$AgreeScore[sum$AgreeCode3==1| sum$AgreeCode2==1| sum$AgreeCode1==1]<-1
table(sum$AgreeScore, exclude=NULL)

#Obtain Prop Agreement.
prop.table(table(sum$AgreeScore)) 

#code 1.
prop.table(table(sum$AgreeCode1))
prop.table(table(sum$AgreeCode2))
prop.table(table(sum$AgreeCode3))

#Cohen's Kappa (factors in chance agreement): 

#For code1 agrement- the way percent agreement is calculated, we get 0.82 agreement.
prop.table(table(sum$AgreeCode1))

#If we calculate it strictly (code1 must match code1x), we get 0.78 agreement.
sum$AgreeCode1<-0
sum$AgreeCode1[sum$Code1==sum$Code1x]<-1
prop.table(table(sum$AgreeCode1))

##calculate Cohen's Kappa.-------
# Contingency table
xtab <- table(sum$Code1, sum$Code1x)
# Descriptive statistics
diagonal.counts <- diag(xtab)
N <- sum(xtab)
row.marginal.props <- rowSums(xtab)/N
col.marginal.props <- colSums(xtab)/N
# Compute kappa (k)
Po <- sum(diagonal.counts)/N
Pe <- sum(row.marginal.props*col.marginal.props)
k <- (Po - Pe)/(1 - Pe)
k #0.7291729 
#0.61-0.80=Substantial
